Al AA-2006-61 44 


Practical Aspects of the 
Equation-Error Method for 
Aircraft Parameter Estimation 


Eugene A. Morelli 

NASA Langley Research Center 
Hampton, VA 


AIAA Atmospheric Flight Mechanics Conference 
August 21-24, 2006 / Keystone, CO 


For permission to copy or to republish, contact the American Institute of Aeronautics and Astronautics, 
1801 Alexander Bell Drive. Suite 500. Reston. VA. 20191-4344 
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NASA Langley Research Center, Hampton, Virginia, 23681-2199, USA 

Various practical aspects of the equation-error approach to aircraft parameter 
estimation were examined. The analysis was based on simulated flight data from an F-16 
nonlinear simulation, with realistic noise sequences added to the computed aircraft 
responses. This approach exposes issues related to the parameter estimation techniques and 
results, because the true parameter values are known for simulation data. The issues studied 
include differentiating noisy time series, maximum likelihood parameter estimation, biases 
in equation-error parameter estimates, accurate computation of estimated parameter error 
bounds, comparisons of equation-error parameter estimates with output-error parameter 
estimates, analyzing data from multiple maneuvers, data collinearity, and frequency-domain 
methods. 
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Nomenclature 

body-axis translational accelerometer measurements, ft/sec 2 
wing span, ft 

mean aerodynamic chord, ft 
covariance matrix 
expectation operator 
mass moments of inertia 

imaginary number = V-T 
cost function 
mass, slugs 

body-axis roll, pitch, and yaw rates, rad/sec 
dynamic pressure, lbf/ft 2 

real part 

standard error 
wing reference area, ft 2 
thrust, lbf 
airspeed, ft/sec 

coordinates of the center of gravity 

coordinates of the reference point 

angle of attack, rad 
sideslip angle, rad 

stabilator, aileron, and rudder deflections, deg 

parameter vector 

variance 
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superscripts 

T 


-1 

f 

subscripts 

o 


transpose 

estimate 

Fourier transform 
time derivative 
matrix inverse 

complex conjugate transpose 


reference value 


I. Introduction 

The equation-error method was one of the first analytical techniques used to estimate aircraft dynamic model 
parameters from flight data 1 . In the 1960’s and 1970’s, the advent of the analog and digital computer, along with 
various theoretical developments and numerous practical successes made the output-error method 2 the 
widely-accepted standard method for aircraft parameter estimation. However, recent analytical developments and 
special problems in aircraft dynamic modeling have highlighted the generality and usefulness of the equation-error 
approach. In this paper, several important practical aspects associated with using the equation-error method for 
aircraft parameter estimation are examined. In some cases, modern treatments of flight data negate earlier problems 
with using the equation-error method for flight data analysis. In other cases, the discussion shows that modern 
implementations of the equation-error method have important advantages for some aircraft modeling problems that 
have become prevalent. 

Following a brief introduction to the equation-error method, the paper is arranged as a sequence of topics, each 
with analytical background and development, and numerical experiments using simulated data. The simulated data 
come mainly from a nonlinear F-16 simulation 3 . Noise sequences similar to those observed in real flight data 
analysis and modeling were added to the simulated aircraft responses. This approach makes it possible to identify 
practical problems and solutions regarding the parameter estimation, because the true parameter values are known 
from the simulation. All of the tools used to generate the results shown in the paper are available in a MATLAB"' 4 
software package called SIDPAC (System IDentification Programs for AirCraft), which is documented in and 
included with Ref. 5. 


II. Aircraft Parameter Estimation using Equation-Error 

In its most common form, the equation-error method calculates aerodynamic parameter estimates that minimize 
the sum of squared differences between values of non-dimensional force and moment coefficients determined from 
measured flight data and corresponding model values. The non-dimensional force and moment coefficients are 
computed by substituting measured and known quantities on the right sides of the following equations 3 
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This produces N values of the non-dimensional force and moment coefficients, where N is the number of data 
points. These values are often called measured force and moment coefficients, even though they not measured 
directly, but rather computed from other measurements and known quantities. 

Model values are computed from a postulated model for each force or moment coefficient in terms of measured 
aircraft states and controls. Although the form of the models can vary, typical linear model structures might be: 
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(10) 


C, = C, f e + c, r jt + c, ±- + C„ t s, + C„ s s , + c„ 


( 11 ) 


For each of the equations (6)-(l 1), substituting force or moment coefficient values computed from Eqs. (l)-(5) 
on the left side, along with corresponding measured states and controls on the right, results in an over-determined set 
of equations for the unknown aerodynamic parameters. This problem can be solved using a standard least-squares 
method. 

For example, the least-squares problem for the pitching moment coefficient C m is formulated using the model 
structure in Eq. (8) as 


Z = XO + v 


( 12 ) 


where 


Z = \C m (l) C m (2) ... C m ( A)] = N x 1 vector of values computed from Eq. (4) 
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= 4x1 vector of unknown parameters 

= TV x 4 matrix of explanatory data vectors or regressors 
T 


v = [v(l) v( 2) ... v(A)l = N x 1 vector of equation 


errors 


The matrix X is assembled using measured data, with each column representing a modeling function, also 
called a regressor. The best estimator of 6 in a least-squares sense comes from minimizing the sum of squared 
differences between the dependent variable measurements z and the model. 
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The least-squares solution for the unknown parameter vector 0 is 

> -1 


0 = 


K4 


X T z 


The estimated parameter covariance matrix is computed from 


Cov{SJmE (0-0)(0-0'f =a 2 [x T X^ ’s [C y ] i,j = 1,2 
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where the number of unknown parameters n p = 4 for this example. The standard errors of the estimated parameters 
are given by the square root of the diagonal elements of the covariance matrix. 


( 4 )* V ^ J ='- 2 ", 


(17) 


There are two important unstated assumptions underlying the expressions given above. One is that the 
regressors, or columns of X , are assumed to be known without error. The second is that the noise on the dependent 
variable z is assumed to be white and Gaussian. In practical aircraft problems, neither of these assumptions holds. 
The consequences of this are discussed below, along with practical ways to correct the problems. 

To apply equation-error parameter estimation in the frequency domain, the data is first transformed from the 
time domain into the frequency domain. For parameter estimation purposes, this can be done very effectively using a 
high-accuracy chirp-z Fourier transform, with the capability to use arbitrary frequencies for the transformation 5 ' 6 . 
Prior to the development of this technique, it was not practical to apply the equation-error method in the frequency 
domain, because the standard method for Fourier transformation of sampled data introduced errors in the 
transformed data and the spacing of the transform frequencies in the frequency band of interest was often coarse. For 
the present work, the frequencies used for all Fourier transformations were 

/ = [0.20 0.24 ... 1.16 1.2] r Hz (18) 


In general, frequency resolution of 0.04 Hz works well. The lower bound should be chosen as 2/T , where T is the 
time length of the maneuver, and the upper bound should be chosen to include the dynamics of interest, typically 
around 1 .2 Hz for rigid-body modes of full-scale aircraft. 

The resulting set of over-determined equations is of the same form given in Eq. (12) above, except that the 
number of data points is m, which is the number of selected frequencies for the Fourier transformation, and the 
problem now involves complex numbers. The analogs of Eqs. ( 14)-( 1 7) for equation-error in the frequency domain 
are 5 
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Note that Eqs. (19)-(22) are the same as Eqs. ( 14)-( 1 7) when the data are real. Consequently, Eqs. (19)-(22) can 
be used for either real or complex data. When the data are complex, the expressions in Eqs. (19)-(22) are equivalent 
to using Eqs. ( 14)-( 1 7) with the real and imaginary parts of the complex data arranged so that the imaginary values 
are stacked below the real parts, and treating the resulting problem as if the data were real. For example, if 
X = X R + jXj , where X R and Xj are N xn p matrices of real numbers, then 


Re(x f x) = Re f (X R - jX, f (X R + jX, ) 


= x$x R +xJx I = 


~X R 

T 


X,_ 


X,_ 


(23) 


and similarly for Re{jc f z\ ■ It follows that the regression problem for complex data can be treated as a problem with 

real data by simply stacking the real and imaginaiy parts of the complex data. However, the value of N in Eq. (16) 
must equal the number of data points, which will be the number of frequency points m . 

The equations above show that the equation-error method for parameter estimation has a relatively simple, 
non-iterative solution, based on linear algebra. The modeling can be done using any of the Eqs. (6)-( 11) (or 
analogous versions of these equations, with different model structures), one at a time. The equation-error method 
can also be applied to state space models to estimate dimensional stability and control derivatives, see Ref. 5. The 
equation-error method can be used very efficiently for very large data sets, e.g, from wind tunnel tests or multiple 
flight test maneuvers. This efficiency is also needed when identifying nonlinear models, which typically involves 
evaluating many different candidate model structures by estimating the parameters for each candidate model 
structure and using statistical modeling metrics to choose among them. The efficiency of the equation-error method 
makes it the only practical choice for this problem. 

The equation-error method can be considered a method wherein the model matches state time-derivative 
information from the dynamic system, rather than matching the states or outputs, as in the output-error method. This 
can be seen from the right sides of Eqs. (l)-(5), where the quantities used as dependent variables in the equation- 
error method are computed from translational and angular accelerations. Consequently, a model that uses equation- 
error parameter estimates will not produce the best match to the measured outputs of the dynamic system, since that 
is not what is being optimized. On the other hand, there is no need to integrate equations of motion to get model 
outputs when using the equation-error method, because the matching is done in the equations of motion themselves 
(hence the name “equation-error”). One important practical consequence is that the equation-error method can be 
applied equally well to data from inherently unstable aircraft flying under closed-loop feedback control. 


III. Practical Issues 


A. Differentiating Noisy Time Series 

Eqs. (3)-(5) show that angular accelerations p, q, r, are required to compute the measured values of the moment 
coefficients C/ , C m , C n . Usually, angular accelerations are either not measured or else the measurements are of poor 
quality. Consequently, angular accelerations must be determined by numerically differentiating the measured 
angular rates, which are noisy. Signal distortion of any kind can degrade modeling results, so it is important that this 
data processing step be done accurately and without introducing additional noise or signal distortion. There are 
many different methods for numerically differentiating noisy data, particularly in the biomedical field 7 . Methods that 
have been found to work well for aircraft data are described here. 
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The main problem with differentiating noisy data is that straightforward finite difference formulas cannot be 
used, because noise on the data causes unacceptable distortions to the relatively small values computed for the 
numerator and denominator in a typical finite difference formula. Consequently, the data processing must be a 
combination of smoothing and differentiation. 

One method is to sequentially fit local approximating polynomial models to the data, then differentiate the 
approximating polynomial. For example, 5 local points (two on either side of the data point where the derivative is 
to be computed) can be fit with a second order polynomial in time, 

y = at 2 +bt + c (24) 


where t is defined as zero at the center point and the constants a, b, and c are estimated using least squares, as in 
Eq. (14). This corresponds to assuming that the acceleration is constant for these five consecutive data points. Then 
the derivative at the center point is computed as the derivative of the local approximating polynomial evaluated at 
the center point where t = 0 , which equals b . This approach can be moved along the time line, and modified for the 
endpoints, where there are fewer than 2 data points on one side of the point in question. References 5 and 8 explain 
this method, and provide practical algorithms for carrying it out. However, in practical cases, the degree of 
smoothing from this approach is often not sufficient, and additional data smoothing must be applied. This can be 
done with zero distortion and zero phase shift using a method based on Fourier analysis and the Wiener filter, called 
an optimal Fourier smoother, which is described in Refs. 5 and 9. The optimal Fourier smoother is a global method, 
because the technique is applied to all of the data from the maneuver at once, rather than operating sequentially on 
local subsets of the data. 

Alternatively, the data can be smoothed initially using the optimal Fourier smoother, then differentiated using 
either the local polynomial method described earlier, or by computing the derivative of the smoothed signal in the 
frequency domain, based on the results from the optimal Fourier smoother. In summary, the three methods are: 

1) Optimal Fourier smoother followed by the local polynomial smoothing derivative 

2) Differentiating the smoothed signal from the optimal Fourier smoother in the frequency domain, then 
reconstructing the smoothed derivative in the time domain 

3) Local polynomial smoothing derivative followed by the optimal Fourier smoother 


To evaluate these methods, random colored noise sequences were used. The objective is to differentiate the 
band-limited portion of the signal, which simulates an arbitrary aircraft response, while rejecting the wide -band 
noise. This situation is very close to practical aircraft 
cases, and covers many possibilities by using random 
noise sequences to generate the band-limited part of 
the signal that simulates the aircraft motion. The 
band-limited part of the signal is known, and its 
derivative can be computed very accurately using 
numerical means. The numerical differentiation 
methods to be evaluated operate only on the noisy 
data, which include both the simulated deterministic 
part and the wide -band noise. 

Figure 1 shows results from applying the three 
methods to 100 different colored noise sequences. 

Methods 1 and 2 gave similar results, with 
root-mean-square errors averaging around 5 percent. 

Method 3 was less accurate, because differentiating 
first (before optimal Fourier smoothing) raised the 
noise level somewhat, and this occasionally led to 
increased errors. By using the optimal Fourier 
smoother first, the noise was mostly removed, and the 



Figure 1. Root-mean-square error for numerical 
differentiation of noisy time series 
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subsequent smooth numerical differentiation was more accurate. Figure 2 shows typical results using all three 
methods on a noisy time series. All three methods approximate the true derivative well, with relatively small 
differences between the true derivative and the numerical derivatives, shown in the plots on the right. 
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Figure 2. Typical example of numerical derivatives of noisy time series 


All of the methods lose accuracy near the endpoints of the time series, but this problem can be eliminated by 
including an extra 1 sec of data on either end and discarding the results for the extra data. This was done for all of 
the results shown in Figs. 1 and 2. The recommendation of 1 sec on either end worked well for 50 FIz data with 
deterministic signal in the range of [0, 1] FIz. These frequencies are typical for practical aircraft problems. The time 
length for an effective endpoint buffer could be different for other cases. 

The optimal Fourier smoother automatically computes a cut-off frequency for the smoother based on the data. 
This cut-off frequency can also be chosen manually by the analyst, based on plots of the Fourier coefficients. The 
best results were obtained when the cut-off frequency was about 0.5 FIz above the frequency band of the 
deterministic signal. This can be easily found visually, based on plots of the Fourier coefficient magnitudes 
produced in the optimal Fourier smoothing method. Again, this recommendation was based on 50 FIz data with a 
deterministic signal bandwidth of [0, 1] FIz. All of the results shown in Figs. 1 and 2 used automatic determination 
of the cut-off frequency, based on the data. 


B. Maximum Likelihood Estimates 

It is sometimes argued that output-error parameter estimates are superior to equation-error parameter estimates 
because output-error parameter estimates maximize a likelihood function, which is a statistical measure of 
conditional probability. Since the output-error parameter estimates have this theoretical underpinning, they are 
considered to be more valid and/or more accurate. In reality, either equation-error or output-error parameter 
estimates can be considered maximum likelihood estimates. The difference between the two approaches is in which 
aspect of the physical system response is being matched by the model, and in the assumptions about the noise. 
Consequently, the true distinction between the output-error and equation-error methods is in which assumptions are 
most useful and valid in modeling the dynamic response of the physical system. 

The basis of maximum likelihood parameter estimation is the likelihood function, which is the conditional 
probability of the measurements given the parameters, assuming Gaussian noise on the measured n 0 outputs 5 , 
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where z,- is the ; th measured output vector, Z represents all the Zj , is the i th model output vector that depends on 
the parameter vector 0 , and R is the measurement noise covariance matrix. For simplicity, assume that there is 
only a single measured output (this is the typical case for the equation-error method), so that Z consists of a vector 
of N measurements taken over the course of a maneuver, and R is a scalar. The maximum likelihood parameter 
estimate is obtained as the value of 6 that maximizes the likelihood function, which is equivalent to minimizing the 
negative log-likelihood function, 


9 = max L(Z; 0) = minl-ln [L(Z; 0)]} = mini — - y t (0)]^ R 1 [z f - y { (0)] + — ln\R\ 1 (26) 

9 o 9 2 . _ ^ 2 I 

where the constant term was dropped because it has no effect on the optimization. For a scalar measurement, 


0 = min \ 


o 2o 7 - 




(27) 


i=i 


where R = cr 2 is the noise variance for the single measurement. 

In the most commonly-used forms, both equation-error parameter estimation and output-error parameter 
estimation are based on minimizing the negative log-likelihood function. The difference is only in the definition of 
the output, and consequently where the noise appears. Assuming the output is an aircraft response, such as pitch 
rate, means that the noise (assumed to be on the model output) appears on the pitch rate measurement. This is the 
output-error approach, which is equivalent to assuming that the equations of motion are deterministic (no process 
noise), with measurement noise on the aircraft response variables. If instead the output is assumed to be 
non-dimensional pitching moment coefficient, the noise (still assumed to be on the model output) appears on the 
non-dimensional pitching moment coefficient. This is the equation-error approach, which is equivalent to assuming 
that the equations of motion include process noise, but there is no measurement noise on aircraft states and controls 
used to form the regressors. Both equation-error and output-error maximize a likelihood function and are therefore 
maximum likelihood parameter estimation techniques. It follows that the theoretical properties 5 of the parameter 
estimates are the same, assuming that the underlying assumptions are valid. The real issue is which set of 
assumptions is most appropriate for the problem, how close reality is to those assumptions, and which aspect of the 
physical system dynamics should be matched by the model, i.e., what should be defined as the measured output to 
be matched by the model. Similar arguments apply for the case of multiple measured outputs. 

A general technique, commonly called filter-error, can be used to find maximum likelihood parameter estimates 
assuming both process noise and measurement noise. However, there are computational and practical problems with 
this approach that have limited its use in practice, see Ref. 5. 

In both the equation-error and output-error methods, the noise covariance matrix R is found using a relaxation 
technique, where 6 and R are optimized alternately while the other is held fixed, and the process is repeated until 
estimates of both 0 and R converge. The equation-error method is normally applied to one force or moment 
coefficient (i.e., one equation) at a time, in which case the relaxation technique is not necessary, because the solution 
from Eq. (27) does not change for a single measured output, regardless of the constant weighting cr 2 . However, the 
estimate ef 2 is necessary for calculation of the estimated parameter covariance matrix, cf. Eq. (16). 

For the single measurement case in Eq. (27) where R = <j 2 , the maximum likelihood estimate of cr 2 is found by 

setting the partial derivative of the negative log-likelihood function with respect to cr" equal to zero, and solving for 

2 

cr , 
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(28) 


i = 1 

Usually, equation-error methods use N — n p rather than N in the denominator above (cf. Eq. (16)), which 
theoretically produces an unbiased estimate of the noise variance, see Ref. 10. For values of N and n p typical of 
flight test maneuvers, N » n p , so this change makes virtually no difference in practice. Similar operations for the 
multiple measurement case give 


( 29 ) 

i=l 

„ 2 

If cr“ or R are not computed from Eq. (28) or (29), then the resulting parameter estimates are not maximum 
likelihood estimates. The practical consequence of this is usually that the estimated parameter standard errors are not 
accurate. 

C. Biased Parameter Estimates 

When the model structure includes either too many terms or too few terms relative to the information content in 
the data, the resulting parameter estimates will be biased. This is true regardless of the parameter estimation method 
employed. In the equation-error method, the fact that the regressors are noisy (because they are assembled from 
measured data) also leads to biases in the parameter estimates. This can be demonstrated by assuming the matrix of 
regressors X is comprised of the true (noise-free) regressors X, , plus a zero-mean random noise vector e , 

X = X, + e (30) 

From Eqs. (12) and (14), the expected value of the estimated parameter vector 0 is computed as 


E 


[0) = E (a 7 *) 1 X T (X,0 + v) 


(31) 


Assuming v and e are uncorrelated with zero mean values, then X = X, + e and v are also uncorrelated, so that 


e[o) = e (a 7 x) ' 


x T x t e 


[d) = E (a 7 a) 1 X t (X-s)0 


E[e) = e-E (* r *) 1 


X T £ 0 


(32) 


The parameter estimate 0 computed using noisy regressors is biased by the quantity -E 




-i 


X T £0 


An F-16 simulation was used to investigate the link between noise on the regressors and bias in the estimated 
parameters using the equation-error method. Geometry and mass properties of the F-16 are given in Table 1. To 
isolate the effects of regressor noise, a linear model for the F-16 short period dynamics was extracted from the 
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nonlinear simulation at a flight condition of 5 deg trim angle of attack and 10,000 ft altitude, using central finite 
differences. The form of the linear model is: 


-Z a a + (\ + Zq)q + Z Ss 8 s 

(33a) 

q = M a cc + M q q + M s S s 

(33b) 
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(33c) 
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This linear model was used to generate simulated data, to avoid any bias in the parameter estimates from model 
structure errors, and thereby isolate the effects of regressor noise on estimated parameter bias. The stabilator input 
was a sequence of square wave pulses with pulse widths in the ratio 2- 1-1-2, where the 1 pulse had a width of 1 sec, 
and the amplitude of the input was 1 deg. Other square waves and frequency sweep inputs with various input 
amplitudes were also used, with similar results to those given next. 

Several methods were used to estimate the non-dimensional model parameters: 

1) Equation-error least squares using noise-free regressors 

2) Equation-error least squares using noisy regressors 

3) Equation-error least squares in the frequency domain using noisy regressors 

4) Equation-error least squares with noisy regressors smoothed using optimal Fourier smoothing 

5) Output-error in the time domain with noisy outputs 

Repeated maneuvers were run using the same input and simulation, but different noise sequences, and all of the 
methods listed above were applied to estimate the model parameters from the simulated data. Figure 3 shows a 
summary of the results. Only the bias is shown for each estimated parameter, which was computed as the difference 
between the mean value of the estimates from the repeated maneuvers, minus the known true value of the parameter. 
This removes the random variation in the parameter estimates. Results were based on estimated parameters from an 
ensemble of 20 repeated maneuvers. 

All added noise was zero-mean Gaussian, with amplitudes typical of flight test data. In general, when the output 
noise levels (i.e., noise levels on C z and C m for the equation-error methods, or on a , q . and a. for output-error) 
were increased, there was an increased variability in the parameter estimates, but no effect on the biases. As 
expected, using least squares with noise-free regressors resulted in unbiased parameter estimates, whereas least 
squares with noisy regressors produced biased parameter estimates. If the analysis was done in the frequency 
domain, or the regressors were smoothed in the time domain first, the unbiased character of the parameter estimates 
was recovered, but the variability in the parameter estimates increased compared to least squares with noise-free 
regressors. The output-error method showed slight biases, but low variability in the parameter estimates for different 
noise sequences. In general, there was good agreement between the results from the output-error method and all of 
the equation-error methods, except for least squares with noisy regressors. Increased noise level on the regressors 
caused increased bias in the estimated parameters from the least squares using noisy regressors, but the other two 
equation-error methods still gave unbiased estimates. Regressor smoothing in the time domain and frequency- 
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domain analysis were almost equally effective, but the frequency domain method showed slightly better 
performance and is computationally more efficient. However, the frequency-domain method does not provide bias 
term estimates, because the data must be detrended prior to applying the Fourier transform. 
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Figure 3. Regressor noise effects on estimated parameter bias 


The frequency-domain method and the smoothing method essentially repair the bias in the estimated parameters 
seen when using least squares with noisy regressors. Both achieve this by smoothing the regressors, since the 
frequency-domain method implements smoothing by transforming the data for a limited frequency band that 
includes frequencies of interest, and the smoothing method does something similar by removing high-frequency 
noise from the time-domain regressors. 


D. Output Matching using Equation-Error Parameter Estimates 

As mentioned earlier, the equation-error method finds parameter estimates that optimize the model match to 
measured derivative data, roughly speaking. However, there is always a desire to know how well a model matches 
measured responses from the aircraft, because they are the most physically meaningful. To make this evaluation 
using equation-error parameter estimates, the meaning of biases must be clearly understood. 

For example, in Eq. (8), if the expansion shown is interpreted as a truncated Taylor series model for the pitching 
moment coefficient, then the modeling functions or regressors should be made from perturbation quantities. If 
instead the directly measured physical quantities are used for the regressors, Eq. (8) becomes 


r =C 

m v ~' n- 


(cr 0 + Her) + C m ^ (<7 0 + Aq ) + C m ^ + A 5 s j + C m ° 


(35) 


where the A quantities are perturbations from the reference values with subscript o. The reference values are 
typically from initial trim conditions. Re-arranging, 
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C m — C m Aoc + C ni Ac/ + C m AS, + 

m m a m q 2 y ^ m§ s s 


Cm a + C'm a a o + C mq Qo + ^m s ^s 0 


(36) 


The bias term in Eq. (36) includes both C m ° and terms related to the steady parts of the regressors, 


C m a„ + 


C mq + Cm,, #s 0 ■ Consequently, the bias term depends on whether the regressors are assembled 


directly from the measured data or from perturbation data where reference values are removed first. The other model 
parameters (stability and control derivatives) are estimated just the same, regardless of whether directly measured 
quantities or perturbation quantities are used. This issue does not occur in cases when the reference values are zero 
for the states and controls used to assemble the regressors, which is common for lateral cases. 


When equation-error parameter estimates are used in the equations of motion to compute aircraft responses, any 
mismatch of the equation-error model to the (derivative) data will be integrated over time. Unless the modeling 
errors from the equation-error parameter estimation are low-amplitude and zero-mean, the result can be a significant 
mismatch between measured aircraft outputs and the model outputs computed from the equations of motion using 
equation-error parameters. In the output-error case, the mismatch between measured outputs and model outputs is 
intentionally minimized by adjusting the parameter estimates. Consequently, aircraft model responses computed 
using equation-error estimates will never show as good a match to the measured outputs as output-error parameter 
estimates. Similarly, using output-error parameter estimates to compute force and moment coefficients will produce 
inferior fits to the data, compared with using equation-error parameter estimates, which are chosen specifically to 
optimize the model fit to force and moment coefficient data. 


Often, analysts will use equation-error estimates to compute model outputs by fixing the equation-error 
parameter estimates and letting the bias terns be estimated using output-error parameter estimation. This essentially 
takes any modeling error that occurs in the equation-error model and assigns it to the bias term. In general, this does 
not work well, and the resulting model fits to measured outputs are poor. However, it is incorrect to indict the 
equation-error parameter estimates for this, because they are being used out of context. It is more meaningful to 
compare the equation-error parameter estimates to output-error parameter estimates, accounting for error bounds and 
any a priori information, and make judgments of the parameter estimate quality based on that information. 


r 

IT. 


A related practical issue occurs when the flight condition drifts slightly from the beginning to the end of the 
maneuver. This affects longitudinal cases where the short-period approximation is being used for the data analysis 
and modeling. The short-period approximation assumes 
a constant airspeed, but often the airspeed is not 
constant, and this can show up roughly as a trend in the 
regressors, such as angle of attack and stabilator 
deflection. The pitch rate can be approximately zero at 
both the beginning and end of a maneuver, whereas the 
angle of attack and stabilator deflection might be 
different at the beginning and end of the maneuver, due 
to flight condition drift and the associated airspeed 
change. The equation-error model for pitching moment 
coefficient will then exhibit a trend over time, as a result 
of including regressors (e.g., angle of attack and 
stabilator deflection) which also have trends over time, 
due to the flight condition drift. The parameter 
estimation algorithm tries to account for the flight 
condition drift in the regressors by adjusting the bias 
term. This causes a visible mismatch with opposite 
signs at the beginning and end of the maneuver, as 
shown in Fig. 4. The solution to this problem is to 
include an airspeed term in the regression model, then 

proceed with the parameter estimation as usual. This is equivalent to removing the flight condition drift, or 
equivalently, applying a time-varying bias corresponding to the airspeed change. In the frequency domain, this 



time, sec 


Figure 4. Pitching moment coefficient model with 
flight condition drift 
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problem is avoided, because any flight condition drift is removed when the data are detrended prior to Fourier 
transformation. 

E. Estimated Parameter Errors 

Because of model simplifications (e.g., linearizing the nonlinear equations of motion) and aerodynamic model 
structure inaccuracies (e.g., approximating the real nonlinear dependencies with tractable mathematical functions of 
instantaneous states and controls), most practical applications of aircraft parameter estimation produce residual 
sequences that resemble colored noise. Typically, a significant portion of the noise power lies in the same frequency 
band as the rigid-body dynamics. A method that corrects for this, giving accurate parameter error bounds for 
time-domain analysis is described in Refs. 5 and 1 1 and outlined here. The method accounts for colored noise by 
modifying the calculation of the covariance matrix in Eqs. (15) and (16). 

The parameter covariance matrix for colored noise is derived by substituting for 0 from Eq. (14) into Eq. (15), 
and using Eq. (12), 


Cov(d) = E ( 9-9)(0-0 ) r = [x T x} * X T E{vv T ^x[x T 


When v is a zero-mean, weakly stationary random process, 

e[vv t ^ = e\vi u, ] = 7Z n , (i~ j) = 7Z VV (j-i) i,j = \,2,...,N 


(37) 


(38) 


where 7Z n , (i — j ) is the autocorrelation matrix for the residuals. The estimated parameter covariance matrix for 
correlated residuals can be computed by substituting for e(vv T j from Eq. (38) into Eq. (37), using an estimate of 
7Z n , (i - j ) computed from 


Kv (*) 


1 

N 


N-k 

I 


i=l 


v i+k 


= K,(-k) 


£ = 0 , 1 , 2 ,.. .,/- 


(39) 


where v ( - is the /* residual. The index £ represents the time separation of the residuals in the summation, and r is the 
maximum time index difference. Since only proximate residuals are significantly correlated, the value of r can be 
relatively small, which reduces the required computations. Combining the last three equations, 


Cov(0) = (a j 'a:) 1 


N N 


y, x i j) x 

. 1=1 7=1 


(* r *r 


(40) 


where xj is they* row of the X matrix, containing the measured regressors at they* data point, and 7Z m {i-j ) is 
replaced by its estimate from Eq. (39). Note that if the residuals are uncorrelated, then 

E(w T ) = ft vv (0)l*a 2 I (41) 


and Eq. (40) reduces to Eq. (15). 

Eqs. (40) with (39) represent a post-processing of the residuals from a least-squares solution to account for the 
effect of correlated residuals on the parameter covariance matrix. The standard errors for the estimated parameters, 
corrected for correlated residuals, are found as the square root of the diagonal elements of the covariance matrix 
calculated from Eq. (40). 
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Figure 5 shows a typical colored noise 
sequence, which is the residual sequence from 
modeling the pitching moment coefficient for a 
doublet sequence using the nonlinear F-16 
simulation with added white Gaussian noise. 

Even in the time-domain plot of Fig. 5, it is 
easy to distinguish prominent low-frequency 
components, which are attributed to 
deterministic modeling error. The approach 
described here assumes that the residual is 
composed of only random components, but the 
low-frequency components are serially 
correlated in time, resulting in colored noise. 

This does not correspond exactly to the 
physical situation, because the low-frequency 
components of the residual are really from 
deterministic modeling errors. Flowever, the 
shift in how the residual is viewed makes it 
possible to compute accurate error bounds for 
the estimated parameters. In effect, the 
problem is changed to the following: given that the true model structure is the one postulated, and the noise is 
correlated in time, what are the errors in the parameters? The fact that the model structure is not exactly correct is 
manifested by colored noise and increased error bounds for the parameters in the specified model structure. 

In the frequency domain, the residuals are no longer serially correlated, because each data point in the regressor 
matrix has changed from measurements at a single sampling time to Fourier coefficients for a single frequency. So, 

the assumption of uncorrelated residuals is 
acceptable in the frequency domain, but the 
assumption that the noise variance magnitude is 
constant over all frequencies (cf. Eq. (21)) is 
not. Figure 6 shows that the power spectral 
density of the model residuals in Fig. 5 does not 
have constant magnitude. However, if the 
frequency-domain analysis is limited to the 
frequency band containing the rigid-body 
frequencies, which is where the large 
components appear in Fig. 6, then the residuals 
can be approximated as having constant 
magnitude. In that case, the theory matches the 
application, and the estimated parameter error 
bounds are calculated accurately without any 
corrections. As a consequence, if equation-error 
is applied in the frequency domain with the 
transform frequencies chosen to approximately 
correspond to the large magnitude components 
of the residuals in the frequency domain, then 
standard calculations outlined above in Eqs. (20)-(22) can be used without modification to calculate accurate 
standard errors for the estimated parameters. 

Computed error bounds are accurate using the frequency-domain method or using the colored noise calculation 
with the regressors smoothed in the time domain. If the regressors are noisy and not smoothed, then using the 
colored noise correction will not give accurate error bounds. 
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Figure 6. Pitching moment residual power spectral density 
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F. Multiple Maneuvers 

Often there is a desire to estimate aerodynamic model parameters using data from more than one flight test 
maneuver. This is useful when individual maneuvers have good information content for estimating some parameters 
but not others. Sometimes multiple maneuvers taken together can provide more complete information for parameter 
estimation. Equation-error allows concatenation of the data from various maneuver by simply stacking the data 
together in the measured output vector z , , and the regressor matrix X in Eq. (12). The implicit assumptions made 
in doing this, in addition to the usual equation-error assumptions, are: 1) the noise v is of the same character for all 
concatenated maneuvers, and 2) the experimental conditions of the concatenated maneuvers are sufficiently similar 
that the model parameters can be expected to be the same. Violating the first assumption leads to inaccuracies in the 
estimated parameter uncertainties. If the second assumption is violated, the parameter estimation results are some 
weighted average of the results that would have been obtained from each of the concatenated maneuvers analyzed 
individually. Parameter estimates from the equation-error method are independent of the ordering of the data points, 
so it does not matter how the data points (i.e., the rows) are ordered in z, and X as far as the parameter estimates 
are concerned; however, the ordering is important for accurately computing parameter standard errors, as discussed 
in the previous section. 

An equivalent method for analyzing multiple maneuvers is to use a priori parameter estimates in the 
equation-error formulation. If the data from two maneuvers are designated by subscripts 1 and 2, and the data from 
both maneuvers are analyzed together, then Eq. (12) becomes 


"*l" 

= 

■*r 

0 + 

"V 

_^2_ 


x 2 _ 


V 2 _ 


(42) 


and the analysis proceeds as usual, assuming that the noise sequences vq and v 2 are similar, and the matrices X l 
and X 2 contain analogous columns of regressors of the same form, implementing a common model structure for 
both maneuvers. If instead the analysis of the data from maneuver 1 was done previously, the cost function to be 
minimized for maneuver 2 can include a priori information about the parameters from the maneuver 1 results, as 
follows 


j{e) = ^(z 2 -x 2 e) T (z 2 -x 2 e)+Ue-0,) T x[\e-e l ) 

2<rf 2 


(43) 


where 


0^ = 


(xfXiY'xfn 


Cov(§i ) = erf (X\ X\ )”' = 


The parameter estimation results in this case are obtained from 

,-i-lr 


9 = 


(x t 2 X 2 )/cj 2 1+ Z? 


(xlz 2 )/a 2 2+ x;'o ] 


Substituting from Eq. (44), 

0 = [(xl X 2 ) I /of + (x[x t ) /of ] [(xl z 2 ) I /of + (x [ Zl ) 1/ 


'2 

a \ 


2 2 2 

Assuming the noise sequences vq and v 2 are similar, erf ~ erf = cr" , so that 


0 = [X T 2 X 2 + xlX l ] [X T 2 z 2 + X[z , " 


(44) 


(45) 


(46) 


( 47 ) 
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which is the same estimate that would be obtained by concatenating the data. Similarly, the covariance matrix 
associated with the parameter estimates in Eq. (45) is 


Cov(e) = [(x$X 2 )/of + 2T 1 ]" 1 = [(*1*2 + [x T , X \ ) 


. n-l 


/CTj 


( 48 ) 


or, assuming the noise sequences V\ and v 2 are similar, 


Cov(t?) = cr 2 T Xl X 2 + X 1 T X 1 


-1 


(49) 


where cr 2 is estimated from Eq. (16), based on all of the data from both maneuvers 1 and 2. Again, this expression 
is the same one that would be obtained by concatenating the data. 

This development shows that as long as the assumptions listed above are acceptable, using a priori values for the 
parameter estimates and covariance matrix is equivalent to concatenating the data from more than one maneuver. 
The result shown here obviously extends to more than two maneuvers. Using a priori values for parameter estimates 
and the covariance matrix is much more convenient than concatenating data, especially for many maneuvers. This 
provides an easy way to incorporate the information from additional maneuvers on parameters already calculated 
from other data. If the noise sequences for the maneuvers to be combined have different variances, the more general 
expressions in Eqs. (45) and (48) can be used. 

In the frequency domain, if multiple maneuvers are transformed using the same frequencies, then the data for the 
combination of multiple maneuvers is simply the addition of the frequency-domain data at each corresponding 
frequency. The analysis then proceeds as usual using Eqs. (19)-(22). This straightforward approach assumes that 
common transform frequencies can be used for multiple maneuvers, which may have different time lengths. This 
can be done with high accuracy using the Fourier transform method described in Ref. 5 or 6. 


G. Data Collinearity 

Data collinearity is a linear or nearly-linear relationship among regressors, which simply means that some model 
terms are nearly the same as a linear combination of other terms in the model. Existence of data collinearity can 
compromise parameter estimation results, because when regressors are nearly collinear, there is ambiguity in how 
the variations in the measured output can be modeled. This is true regardless of the parameter estimation method, 
but the issue is easier to study in detail for the equation-error method. 

For the equation-error method, the detrimental effects of data collinearity are closely tied to the noise levels on 
the measurements. In fact, if the noise levels are zero, data collinearity has no effect at all on least-squares parameter 
estimates. This was demonstrated using simulated data for longitudinal short period linear dynamics of the F-16. 

A doublet sequence was applied to the stabilator, and the simulated aircraft responses were corrupted with white 
Gaussian noise. A simple pitch rate feedback to the stabilator was used to correlate the stabilator and pitch rate 
regressors. This situation is a common cause of data collinearity in practice. By adjusting the feedback gain, the 
regressors could be correlated to an arbitrary extent. Pair-wise regressor correlations were defined by: 


C 


r jk 


jk 


fCfk 


Cox 

\0jA\ 

) 

yj Var (i 

°j) 

Var ( 

' 4 ) 


j,k = \,2,...,n p 


(50) 


-\<r jk <\ 


j,k = 1 , 2 ,...,« 


( 51 ) 


With no noise on the measured output ( C z or C m ) and using noise-free regressors, the ordinary least squares 
solution given in Eqs. ( 1 4)-( 1 6) very accurately computed the true parameter values, even when the maximum 
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absolute value of the regressor correlation exceeded 0.99. The solid line in Fig. 7 shows this result for the C m 
parameter. 
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This demonstrates that extremely high data collinearity can be handled using the default double-precision 
arithmetic in MATLAB®, as long as there is no noise. The problems occur when measurement noise is added. 
Keeping the regressors noise-free and adding noise to the measured output produces approximately unbiased 
estimates, as shown in Fig. 7. Flowever, when 
noise is added to the regressors, the estimated 
parameter bias is significant. If least squares in 
the frequency domain is applied to the case 
where both the regressors and the measured 
output are noisy, nearly unbiased parameter 
estimates are recovered, and the parameter 
estimates are as good as when the regressors had 
no noise at all. This result comes from using a 
limited frequency band for the Fourier transform, 
as discussed earlier, which effectively removes 
the wide -band noise from the regressors in the 
frequency domain. Similar results were found for 
the other model parameters, and for inputs 
different than a doublet. 
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In Fig. 7, note that when the maximum 
regressor correlation exceeds approximately 0.9, 
the results degrade in all cases. Consequently, a 
good practical rule of thumb is that the effects of 
data collinearity get worse when the absolute 
value of regressor correlation exceeds 0.9. Increased noise levels show the same trends, but the bias in the estimated 
parameters increases for least squares with noisy regressors. Increased noise on the outputs does not affect the bias, 
but increases the random variation in the parameter estimates in all cases. This empirical finding is consistent with 
the expression derived for the bias in Eq. (32). 


Figure 7. Data collinearity and noise effects on estimated 
parameter bias 


In summary, noise on the regressors combined with data collinearity will produce significant biases in estimated 
parameters using least squares in the time domain. Measurement noise must be present as a catalyst for the data 
collinearity to cause inaccuracies in the parameter estimates. If the regressors and outputs are noise-free, least 
squares using ordinary double-precision arithmetic can accurately estimate model parameters for pair-wise regressor 
correlations as high as 0.99; however, when noise is added to the regressors, the estimated parameters are biased. 
Using least squares in the frequency domain repairs the adverse effects of data collinearity by effectively removing 
the wide -band noise on the regressors. 


IV. Conclusions 

Several practical issues related to applying the equation-error method to aircraft parameter estimation problems 
were examined and discussed. The issues included numerical differentiation of noisy data, maximum likelihood 
parameter estimation, biases in estimated parameters, computing aircraft responses using parameters estimated with 
equation-error methods, calculating estimated parameter errors, parameter estimation based on data from more than 
one maneuver, and data collinearity. Some limitations and deficiencies of the equation-error method were shown to 
be removed by advances in applying the technique. Specifically, applying the equation-error method in the 
frequency domain was shown to have advantages in computational efficiency, accurate computation of the estimated 
parameter errors, minimizing biases in the estimated parameters, and robustness to the adverse effects of data 
collinearity. The equation-error method is applicable to aircraft modeling problems that are increasingly prevalent in 
modern aircraft, such as nonlinear aerodynamic modeling, large data sets, and aircraft that are open-loop unstable. 

Practical issues were examined using data from an F-16 nonlinear simulation, with realistic colored noise 
sequences added to the computed aircraft responses. This allowed a clear view of the effectiveness of various 
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approaches to modeling problems, because the true values of the model parameters were known. This approach is 
necessary to fully understand and validate tools for flight data analysis and modeling. 

Software tools for all the techniques and investigations included here, along with many other methods used in 
aircraft system identification, can be found in the software package called SIDPAC (System IDentification Programs 
for AirCraft), which is available with Ref. 5. 

The venerable equation-error method, implemented using modem analytical approaches and numerical 
techniques, was shown to be accurate and useful for aircraft parameter estimation problems. 
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Table 1 . Geometry and mass properties 
of the F-16 aircraft 


length c , ft 

11.32 

wing span b , ft 

30 

wing area S , ft 2 

300 

X re f , ft 

0.35 c 

y r ef ’ ft 

0.000 

z ref ' ft 

0.000 

x cg , ft 

0.25 c 

y eg ’ ^ 

0.0 

z cg ’ ft 

0.0 

m , slugs 

647.2 

I x , slugs-ft 2 

9,496 

/ , slugs-ft 2 

55,814 

I. , slugs-ft 2 

63,100 

I xz , slugs-ft 2 

982 
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